home *** CD-ROM | disk | FTP | other *** search
/ Super PC 34 / Super PC 34 (Shareware).iso / spc / UTIL / DJGPP2 / V2 / DJTST200.ZIP / tests / libc / ansi / math / elefunt / tasin.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-06-01  |  6.9 KB  |  313 lines
  1. /* -*-C-*- tasin.c */
  2.  
  3. #include "elefunt.h"
  4.  
  5. /*
  6. #     program to test asin/acos
  7. #
  8. #     data required
  9. #
  10. #        none
  11. #
  12. #     subprograms required from this package
  13. #
  14. #        machar - an environmental inquiry program providing
  15. #                 information on the floating-point arithmetic
  16. #                 system.  note that the call to machar can
  17. #                 be deleted provided the following four
  18. #                 parameters are assigned the values indicated
  19. #
  20. #                 ibeta  - the radix of the floating-point system
  21. #                 it     - the number of base-ibeta digits in the
  22. #                          significand of a floating-point number
  23. #                 irnd   - 0 if floating-point addition chops,
  24. #                          1 if floating-point addition rounds
  25. #                 minexp - the largest in magnitude negative
  26. #                          integer such that float(ibeta)**minexp
  27. #                          is a positive floating-point number
  28. #
  29. #        ran(k) - a function subprogram returning random real
  30. #                 numbers uniformly distributed over (0,1)
  31. #
  32. #
  33. #     standard fortran subprograms required
  34. #
  35. #         abs, acos, alog, alog10, amax1, asin, float, int, sqrt
  36. #
  37. #
  38. #     latest revision - december 6, 1979
  39. #
  40. #     author - w. j. cody
  41. #              argonne national laboratory
  42. #
  43. #*/
  44.  
  45. /* #undef ABS
  46. #define ABS fabs */
  47.  
  48.  
  49. void
  50. tasin()
  51. {
  52.     int i,
  53.         ibeta,
  54.         iexp,
  55.     irnd,
  56.         it,
  57.     i1,
  58.         j,
  59.         k,
  60.         k1,
  61.         k2,
  62.         k3,
  63.         l,
  64.         m,
  65.         machep,
  66.         maxexp,
  67.         minexp,
  68.         n,
  69.         negep,
  70.         ngrd;
  71.     float a,
  72.         ait,
  73.     albeta,
  74.         b,
  75.     beta,
  76.         betap,
  77.         c1,
  78.         c2,
  79.         del,
  80.         eps,
  81.         epsneg,
  82.         half,
  83.         r6,
  84.         r7,
  85.         s,
  86.         sum,
  87.         w,
  88.         x,
  89.         xl,
  90.         xm,
  91.     xmax,
  92.         xmin,
  93.     xn,
  94.         x1,
  95.         y,
  96.         ysq,
  97.         z,
  98.         zz;
  99.  
  100.     machar(&ibeta, &it, &irnd, &ngrd, &machep, &negep, &iexp, &minexp,
  101.     &maxexp, &eps, &epsneg, &xmin, &xmax);
  102.     beta = (float) ibeta;
  103.     albeta = ALOG(beta);
  104.     half = 0.5e0L;
  105.     ait = (float) it;
  106.     k = (int) (ALOG10(ipow(beta, it))) + 1;
  107.     if (ibeta != 10)
  108.     {
  109.     c1 = 201.0e0L / 128.0e0L;
  110.     c2 = 4.8382679489661923132e-4L;
  111.     }
  112.     else
  113.     {
  114.     c1 = 1.57e0L;
  115.     c2 = 7.9632679489661923132e-4L;
  116.     }
  117.     a = -0.125e0;
  118.     b = -a;
  119.     n = 2000;
  120.     xn = (float) n;
  121.     i1 = 0;
  122.     l = -1;
  123.  
  124.     /* random argument accuracy tests */
  125.  
  126.     for (j = 1; j <= 5; ++j)
  127.     {
  128.     k1 = 0;
  129.     k3 = 0;
  130.     l = -l;
  131.     x1 = ZERO;
  132.     r6 = ZERO;
  133.     r7 = ZERO;
  134.     del = (b - a) / xn;
  135.     xl = a;
  136.  
  137.     for (i = 1; i <= n; ++i)
  138.     {
  139.         x = del * ran(i1) + xl;
  140.         if (j <= 2)
  141.         {
  142.         y = x;
  143.         ysq = y * y;
  144.         }
  145.         else
  146.         {
  147.         ysq = half - half * ABS(x);
  148.         x = (half - (ysq + ysq)) + half;
  149.         if (j == 5)
  150.             x = -x;
  151.         y = sqrt(ysq);
  152.         y = y + y;
  153.         }
  154.         sum = ZERO;
  155.         xm = (float) (k + k + 1);
  156.         if (l > 0)
  157.         z = /* (double) --buers */ asin(x);
  158.         if (l < 0)
  159.         z = /* (double) --buers */ acos(x);
  160.  
  161.         for (m = 1; m <= k; ++m)
  162.         {
  163.         sum = ysq * (sum + 1.0e0L / xm);
  164.         xm = xm - 2.0e0L;
  165.         sum = sum * (xm / (xm + 1.0e0L));
  166.         /* xm -= 2.0L;
  167.         sum *= xm / (xm + 1.0L); */
  168.         }
  169.  
  170.         sum = sum * y;
  171.         /* sum *= y; */
  172.         if (!((j != 1) && (j != 4)))
  173.         {
  174.         zz = y + sum;
  175.         sum = (y - zz) + sum;
  176.         /* sum += (y - zz);      */
  177.         if (irnd != 1)
  178.             zz = zz + (sum + sum);
  179.             /* zz += sum + sum;       */
  180.         }
  181.         else
  182.         {
  183.         s = c1 + c2;
  184.         sum = ((c1 - s) + c2) - sum;
  185.         zz = s + sum;
  186.         sum = ((s - zz) + sum) - y;
  187.         s = zz;
  188.         zz = s + sum;
  189.         sum = (s - zz) + sum;
  190.         if (irnd != 1)
  191.             zz = zz + (sum + sum);
  192.         }
  193.         w = 1.0e0L;
  194.         if (z != ZERO)
  195.         w = (z - zz) / z;
  196.         if (w > ZERO)
  197.         k1 = k1 + 1;
  198.         if (w < ZERO)
  199.         k3 = k3 + 1;
  200.         w = ABS(w);
  201.         if (w > r6)
  202.         {
  203.         r6 = w;
  204.         x1 = x;
  205.         }
  206.         r7 = r7 + w * w;
  207.         xl = xl + del;
  208.     }
  209.  
  210.     k2 = n - k3 - k1;
  211.     r7 = sqrt(r7 / xn);
  212.     if (l < 0)
  213.     {
  214.         printf("\fTEST OF ACOS(X) VS TAYLOR SERIES\n\n");
  215.         printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n", n);
  216.         printf("      (" F15P4E "," F15P4E ")\n\n\n", a, b);
  217.         printf(" ACOS(X) WAS LARGER%6d TIMES,\n", k1);
  218.         printf("             AGREED%6d TIMES, AND\n", k2);
  219.         printf("        WAS SMALLER%6d TIMES.\n\n", k3);
  220.     }
  221.     else
  222.     {
  223.         printf("\fTEST OF ASIN(X) VS TAYLOR SERIES\n\n");
  224.         printf("%7d RANDOM ARGUMENTS WERE TESTED FROM THE INTERVAL\n", n);
  225.         printf("      (" F15P4E "," F15P4E ")\n\n\n", a, b);
  226.         printf(" ASIN(X) WAS LARGER%6d TIMES,\n", k1);
  227.         printf("             AGREED%6d TIMES, AND\n", k2);
  228.         printf("        WAS SMALLER%6d TIMES.\n\n", k3);
  229.     }
  230.     printf(
  231. " THERE ARE%4d BASE%4d SIGNIFICANT DIGITS IN A FLOATING-POINT NUMBER\n\n",
  232.         it, ibeta);
  233.     w = -999.0e0;
  234.     if (r6 != ZERO)
  235.         w = ALOG(ABS(r6)) / albeta;
  236.     printf(" THE MAXIMUM RELATIVE ERROR OF" F15P4E " = %4d **" F7P2F "\n",
  237.         r6, ibeta, w);
  238.     printf("    OCCURRED FOR X =" F17P6E "\n", x1);
  239.     w = AMAX1(ait + w, ZERO);
  240.     printf(
  241.         " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS" F7P2F "\n\n\n",
  242.         ibeta, w);
  243.     w = -999.0e0;
  244.     if (r7 != ZERO)
  245.         w = ALOG(ABS(r7)) / albeta;
  246.     printf(" THE ROOT MEAN SQUARE RELATIVE ERROR WAS" F15P4E " = %4d **" F7P2F "\n",
  247.         r7, ibeta, w);
  248.     w = AMAX1(ait + w, ZERO);
  249.     printf(
  250.         " THE ESTIMATED LOSS OF BASE%4d SIGNIFICANT DIGITS IS" F7P2F "\n\n\n",
  251.         ibeta, w);
  252.     if (j == 2)
  253.     {
  254.         a = 0.75e0L;
  255.         b = 1.0e0L;
  256.     }
  257.     if (j == 4)
  258.     {
  259.         b = -a;
  260.         a = -1.0e0L;
  261.         c1 = c1 + c1;
  262.         c2 = c2 + c2;
  263.         l = -l;
  264.     }
  265.     }
  266.  
  267.     /* special tests */
  268.  
  269.     printf("\fSPECIAL TESTS\n\n");
  270.     printf(" THE IDENTITY  ASIN(-X) = -ASIN(X)  WILL BE TESTED.\n\n");
  271.     printf("        X         F(X) + F(-X)\n");
  272.  
  273.     for (i = 1; i <= 5; ++i)
  274.     {
  275.     x = ran(i1) * a;
  276.     z = asin(x) + asin(-x);
  277.     printf(F15P7E F15P7E "\n\n", x, z);
  278.     }
  279.  
  280.     printf(" THE IDENTITY ASIN(X) = X , X SMALL, WILL BE TESTED.\n\n");
  281.     printf("        X         X - F(X)\n");
  282.     betap = ipow(beta, it);
  283.     x = ran(i1) / betap;
  284.  
  285.     for (i = 1; i <= 5; ++i)
  286.     {
  287.     z = x - asin(x);
  288.     printf(F15P7E F15P7E "\n\n", x, z);
  289.     x = x / beta;
  290.     }
  291.  
  292.     printf(" TEST OF UNDERFLOW FOR VERY SMALL ARGUMENT.\n\n");
  293.     x = pow(beta, (float) minexp * 0.75e0L);
  294.     y = asin(x);
  295.     printf("       ASIN(" F13P6E ") =" F13P6E "\n\n", x, y);
  296.  
  297.     /* test of error returns */
  298.  
  299.     printf("\fTEST OF ERROR RETURNS\n\n\n");
  300.  
  301.     x = 1.2e0L;
  302.     printf(" ASIN WILL BE CALLED WITH THE ARGUMENT" F15P4E "\n",x);
  303.     printf(" THIS SHOULD TRIGGER AN ERROR MESSAGE\n\n\n");
  304.     fflush(stdout);
  305.     errno = 0;
  306.     y = asin(x);
  307.     if (errno)
  308.     perror("asin()");
  309.     printf(" ASIN RETURNED THE VALUE" F15P4E "\n\n\n", y);
  310.  
  311.     printf(" THIS CONCLUDES THE TESTS\n");
  312. }
  313.